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NUMERICAL METHODS OF SOLVING A SYSTEM OF MULTI- 
DIMENSIONAL NONLINEAR EQUATIONS OF THE 
DIFFUSION TYPE 


A. V. Agapov, B. I. Kolosov 

it 

INTRODUCTION /3 

The non-steady state' equations of the diffusion type are the most widely 
used equations, to which the description of many important physical phenomena are 
reduced. The literature [5, 6] has made an extensive study of the linear class of 

t 

these equations, and many numerical methods for solving them on computers in the 
one-dimensional and the multidimensional geometry have been developed and have 
become classical. The situation is different in the case of the nonlinear depen- 
dence of the diffusion coefficients on the initial solution, which frequently 
arise when studying nonlinear effects in. the theory of quasi-linear plasma oscil- 
lations, in the problems of heat and mass transfer, and other present day physical 
applications. This article proposes a universal numerical algorithm for solving 
these problems, having the property of conservation and numerical stability. The 
formulation and study of the algorithms presented below are based on a priori in- 
formation connected with the initial problem. 

We -should also note that since certain existing physical phenomena are ex- 
amined, which are modeled by nonlinear equations of the diffusion type, the exis- 
tence for solutions of these equations is a natural requirement, which indicates 
the useful information of the initial physical model. It is known that diffusion 
equations are based on a balance of .th^fjnw of physical quantities determined by 
a differential eiJ^ression of the form , where f is the physical quantity 

examined; — diffusion coefficient. Since in modem physics the laws 
of conservation are only - considered in geometric space, but in the phase space of _/4 

I 

velocity or the Fourier transformation space, in the following .the variable 

, . . . . j . 

. 

Numbers in the margin indicate pagination of original foreign text. 


1 



will designate the parameters of the corresponding space [for example, 
or ] j in. which the diffusion processes are studied. As a rule, 

nvnnerical modeling of non-steady state phenomena is always reduced to discretiza- 
tion of the tame variable for a sequence of evolutionary layers, for which the 
"conservation" of the solution algorithm.is of great importance. This solution 
is the consequence of the necessary agreement hetr^een the a priori requirements 
of the differential problem and the properties of its difference analog. Here, 
"conservation" is reached when, when the stable algorithms are foamiulated, the 
values of f and are taken on the same evolutionary layer, which leads in 

its turn to implicit ' difference analogs of the initial equations with respect to 
the solution and, as a result, to systems of nonlinear equations which may be 
solved only using iteration approaches. 

In conclusion, we should note that since numerical methods of solving multi- 
dimensional equations are. examined here, whose realization falls within the frame- 
work of modem computers, one of the basic requirements imposed on the algorithms 
is a minimum number of computational operations for a given accuracy of'-the ap- 
proximation to the solution. 


1. FOEMULATION AND 'BASIC PROPERTIES OE THE PROBLEM. 


Let us now formulate the problem. For this purpose we shall consider 
Euclidean -space 

\ ' R.. ^ ■ ft 

Let us now assume that QC is a certain region in R^ and let us set' (3,=^ 
and AV ~ r’X'(o^T) is the lateral surface Q^, Now let us examine on a certain 

real function f , having the first and second generalized- derivative. For each 
ie , let us introduce the Banach space and the Sobolev space » 

which consists of the functions of the space t.Q) , having in Q the generalized 

derivatives summed with the square, Wj(Q) produces Hilbert space with respect 

to the norm 


/5 




CD 


Let us now assume f belongs to the class of these functions F in , 

- "T^ 

which for each.-^e£OpTj produces Hilbert space HCQ) either with respect to the 


2 



norm 


w'(a) 

or with respect to its equivalent. 


Under these conditions, let us examine the quasi-linear equation in the 
domain Qt with the divergent right side 


• ^ ^ .2. 1. P . 


( 2 ) 


I. 


for which t is given on the surface S^. Then for the known diffusion coeffi- 
cients the equation (2) determines the unique solution. We shall also assxmie 

that the operator 

7 • ^ 75 — 

4 = 


. A .y 

exists in Hilbert space Li((X) , and the set is its domain of definition, 
consisting of the functions (Q) such that it is degener- 
ate, i.e., for any non-zero vector f inequality is satisfied 




(3) 


where are positive numbers. The condition (3) will be called the condition 

of strong ellipticity. If 


then the condition of strong ellipticity has the form 

where t‘y^ are certain continuous positive functions. 


(4) 16 


Now to complete the formulation of the problem, we must add additional con- 
ditions on equation (2) which determine the dependence of the coefficients 
on the solution f . For this purpose, let us write the following equation in the 
domain 



(5) 
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where in the general case is a certain linear operator which acts in Hilbert 

space H, produced by at give a sufficient criterion 

for ellipticity of the operator / and as a result the parabolic nature of the 

•*( 

system of equations (2) - (5) . Let us assume L exists and is positive. Then 
we have the follox^ing: 

LeTnma 1 . If for tbe inequalities and 

0 ) 

hold, then the matrix is positive definite. 

Actually, it follows from the positive nature of A ' that f 


where K is the cone for Q 




Since according to the conditions of the Lemma 

in a similar way,^^ 


r 



from which we have 


or in view of /pSV » which designates 


the positive definiteness [^<^1 • Although the properties of the smoothness of 

the function , with respect to which we shall assume that the condi- 

tion of the Lipschitz continuity is satisfied at the norm W* (jQ) , mainly 




( 6 ) 


We should note that equation (5) means that from the functional equation determin- 
ing the diffusion coefficients, we may distinguish a certain linear part for . 
Thus, the system of equations (2) and (5) under the ellipticity conditions (3) and 
given on the surface of the boundary conditions completely determines the solu- 


tion and 


n 


Let us turn to a numerical solution of the problem. For this purpose, we 
write the semi-difference analog of the system of equations (2) - (5) , which has 
the following form in view of the conservation requirements; 


f’ ^ ■ 


(7) 


C8) 


where n designates the number of the evolutionary layer and the time layer in 


4 



the given case, and At — the discretization step of the variable &\ 


i T ^ 


Let us examine the elliptical operator in the domain Q 

A 


s - -I SL Tj" 2. r -i 

i' 


C9) 


where is determined on ^-s-T V for . 


Since the solution of the system of equations (7) for each value of the 
evolutionary layer n is reduced to inversion of the elliptical operators g , 
which may be satisfied in view of the nonlinear nature only using iteration ap- 
proaches, we shall now examine several properties of the operators which are not 
omitted ‘(9) . 

' ^ ' 

Let us assume ^ are continuous at P for , and 

^ > 

we shall assinne is a factor of all the finite functions in Hilbert space H, 

which produce 2 £ at • We shall use £ to designate the 

closures Then Is called the minimtini operator producing , The defini- 

tion domain is the space at 0 and ■‘^'1 sr-O . By the expan- 

sion of the operator , we designate the operator at the set 

’.W * )j where ( Q ^ is the set ,^(^J such that 

of the equation 

. I/I- 

is called, as is known, the solution of the first boundary value problem, For it 
the following estimate holds 


X.ir 


*0 


Then, the solution 


'll 


Cillfll,; 


iQ)- 




CIO) 


2. ITERATION SOLUTION OF A SYSTEM OF NONLINEAR ELLIPTICAL EQUATIONS. 

Iteration methods are widely used for inversion of linear elliptical syf3'».' " 
stems of the form (7) - (8) , However, here the requirement for them is dictated 
also by the nonlinearity of the system, and in contrast to the linear case, a 
study of the convergence of the iteration process is of great importance. For 
convenience, we shall examine the initial proximation individually 
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( 11 ) 



( 12 ) 


and the iteration continuation Ox the system (11) - (12) for the solution 

vl~- ( 8 ); 



(13) 

(14) 


The advantage of distinguishing the initial approximation is due to the 
fact that for certain functions , the unbounded operators > and the 

presence of a large amount of a priori information regarding the solution is s^affi- 
cient to confine ourselves to a one’-step process (11) - (12), if the solution of 
the equation (11) is fotmd with a high order of accuracy. This approach, which is 
related to the Introduction of non-control grids, will be examined below. In the 
case of bounded operators , a system of the predictor-corrector type (13) - (14) 
is preferable. We shall study its convergence. 


Let us assume _Li_ is positive definite in H and a bounded operator, i.e., 

. Then the following is valid; 

Theorem 1 . The iteration process (13) - (14) converges for each n to 
the solution of the system (8), if the inverse operator £* is positive and 

the inequalities ^ . 

hold. ' 


/9 


In view of the positive nature of L , there is an inverse operator, and 
/I /T*|j > f3:^om which the solution of the equation (8) follows. 


Let us introduce the notation 






(15) 


Then, subtracting equation (13) from equation C7) , with allowance for the notation 
in (15), we obtain 
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( 16 ) 


> '4 2 A ■ 2- i ” ‘ 

^ T fw) » 


Now let us determine the rigid expansion determined hy the operator 


/'"P*’ _ J.P” 9 *" - 

^ T(^> ■ V . 


9*^. 


Since the matrices 2)^^ and , in view of.the conditions imposed on the 
function , have positive minimum eigenvalues, the inequality (10) follows. 
We have the following from Lemma 1 


and (13) ^ • Consequently ' is bounded. 


Let us also assume 
we obtain 




.y 


Determining (16) 


ii"" [ li&rSwJ t ijfOwi ^ !l AulUl 1. 


In its turn, from (17) we have 

C,(f )llU„tlwi‘ . 


Now, subtracting equation (8) from (14) and applying the condition of the Lipschitz /lO 

, we find 


continuity (6) for the function, ^ 


L = Cjj ( f fm-o , ^ , Xp'*) - ( 17 ^ ) 


(18) 


from which, by stimming the ineqvialities (18) with respect to a and 3 ? we obtain 
the conditions of the theorem, which completes the proof. 


A necessary condition for the implementation of the iteration process (13) 

- (14) is, as is known, its stability with respect to small perturbations. Since 
the stability is due to the correct nature of the difference problem, we shall ex- 
amine a certain general approach to the control of the system (13) - (14) for an 
arbitrary positive-definite and bounded operator L, and then its specific modifica- 
tions as applied to the integral operator L and the differential (tinbounded) opera- 
tor 
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3. GENERAL (CONSERVATIVE) ALGORITHM OF THE ITERATION CONTROL. 


Let us consider the numerical implementation of the process (13) - (14) , 
when L is a positive definite operator. For this purpose, we use the following 
control iteration process. We introduce a certain positive self-conjugate opera- 
tor K in Hilbert space E, in which L exists. Let us assxime K: 


H ❖ H : H, ^5^ H ; Wt : <>® 

, - - 

Let us set Ma= Su\3 (lLx,x') . Let us examine the following iteration process [1] 

iwi^l 


s - . a;.r- 




(19) 

( 20 ) 


If V / belongs to the region of values 4 ? then the following 

theorem is valid [1]. 

Theorem 2 . For ^ the sequence / from (19) - (20) 

converges to the solution of the equation (14) . 

To prove this statement, instead of the space examined H, i.e., 
with ^ , a certain nex<r space is introduced for 

namely, readily shown-that the operator 

is self-conjugate in « Let us set ■ 

and then 


/II 




C21) 


We should note that the norms produced by the scalar products 

are equivalent. It follows from (12) that . The proof of Theorem 2 

follows from this fact and the Krasnosel ’skiy theorem. 


However, Theorem 2 does not make it possible to applybthe process (19) - 
(20) to ■ the solution of the system (13) - (14), In its turn, the possibility of 
a numerical approximation of the process (19) - (20) establishes Theorem 2 IlJ. 
Let us examine the process 


i 





( 22 ) 
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i.e., we assiime the operator LL is known x^ith a certain non-self-conjugate 
error , and we assume the same hold for - We shall assume that 

I > 1 f II 1 f i ^ ® , II ^ 4^ 11 ^ equation 

T(iV/cTdil! <T - 


Theorem 3 , Let 4 ;; he the solution of equation (8) . Then the follow- 
ing estimate holds: 


where 


1 - If !.<«;■ 

0 & 


, then 


'? ■■ -r^. 


2. If q = 1, then 




The proof of Theorem 3 follows from the proof of the preceding estimate 
Let us consider the convergence rate of .the iteration processes (11) - (14) . 


/12 


Here, the following Theorem 3 is valid [2j. If S*”-“ 

»^rs Ci 




belongs to the 


region of values for the operator CL and =0 , then 

!i4?-4;«iKT 


The proof of the theorem is based on the follox^ing 

Theorem 3 may be refined for finite-dimensional space, namely, for finite-dimen- 
sional space as applied to the process (12) . It is known that J are two matrices 
and , so that 

(CL-f-ALyZ-Cai fV4A„ 

Here we may assume that the following estimate is satisfied [ij 
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Then 




where 


■ i*s 


i. 


tntn 

J? - ■' 


Let us refine the form of the matrices Z and Z . . Z is the matrix whose 

' .y.' 1 

columns are K-normed vectors ■ , where Uj ate the eigenvectors of the 

symmetric matrix jS;"* ? P — dimensionality of the matrix; K — the quantity 


A = 0. 


Thus, the solution, of the system (7) - (8) may he reduced to two iteration 
processes: the process (13) - (14) and the process (19) - (20). Both of these 

processes cannot he implemented in practice, since the two-dimensional nature of 
the parabolic equation makes the time for the solution of (7) - (8) beyond the 
limits of present day computers. This article proves the convergence of both 
processes to one, namely, the following processes examined 


■ ' 2. 
^ 


|. 


iw-r 

Z A + K ^ A ( f ) 


s. h'*" • I" = 




(23) 

C24) 

(25) 


4, CERTAIN SUFFICIENT CONDITIONS FOR THE CONVERGENCE OF THE CONTROL PROCESS 
C23) - C25). 

Let us consider the process (23) - (25) , where the ooperator K is selected 
as follows: 


1 . 

2. K ^ acts in the cone, i.e,, it is positive. 



Letvus examine the new operator SK , where V is a given number. We shall 

select £ so that the conditions for the convergence of the process (23) - (25) are 
satisfied with respect to , i.e., this operator changes a cone into a 

cone. In addition, K & K • Under these conditions in space H, X'irhich produces 
W^(Qy , the following is valid: 


Theorem 4 . In the space H examined, the process (23) - (25) converges if 
Let us subtract equation .'(3) from equation (12) 


< 

>Nf ^ - 

It may be readily shoxjn that ' 


. lp,y 


(26) 


satisfies the ellipticity condition 


Let us set 
and 
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, which is readily proven by induction. Actually 


' 7's ° = 


These equations are strictly larger than zero on S(f) . A similar statement is 
valid for 

3 tut 3;r- »• 

since for the first iteration indicated above, the induction is valid. At the same 
time we should note that the minimum eigenvalue of the matrix 

is > Pi 9 > 0^. 

Let us introduce the expansion of .Ljs for the operator 
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1^ 9 - 'i 1 9 


1 ? '__ 7 / ‘ « 

for which, as was shown above, the following estimate is valied (10) 




In addition, from this equation, taking into account (26) , we obtain 




Now, taking the fact into account that the Sobolev spaces or form a 

Hilbert scale of space, we have the inequality 


Let us subtract (.4) from Cl3) . We thus have 

K ( sTf!- 


^ A^.,. f, 




or 


A 




-UJ, 


C27) 


Let us deteinnine the norm (27), from which x^e use the Lipschitz continuity condi- 
tion ' • ‘ . 

where we set 4^ • 
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Finally, combining the inequalities obtained, we have 




xjhich must be proven. 


We should note that everywhere for ■ , the following norm is selected; 

For this norm, all of the computations examined above may be satisfied most simply. 
The proven theorem establishes the convergence of the iteration process in Sobolev 
space. However, certain estimates, for the purpose of distinguishing the control 


12 



process from the non-control process, may be obtained in the metric C, which will 
be given below. 


5 . STABLE ITERATION ALGORITHM FOR SOLVING THE SYSTEM' (13) - (14) , WHERE L IS 
THE INTEGRAL OPERATOR OE THE FIRST KIND, 


Let us examine the system (13) - (14) for a specific form of the operator 
L, namely: 


where Xe (i* V is e certain positive kernel, Then the system Cl3) - 

(14) assumes the form 


where 


•o” P*'*' ’ 


er®’)- 


(28) 

(29) 


The presence of the Fredholm operator of the first kind does not make it 
possible to immediately solve the equations (13) ~ (14) , Therefore, we shall use 
the algorithm (23) - (25) . Let us reduce (29) to a system of integral relation- 
ships, namely, let us introduce in the domain Q the difference nodes and we shall 
designate them by letters in the text, Then for Xj Xj, Q 

we obtain -the following integral relationships; 
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where 



(30) 


The study [2] was devoted to solving the system of the form (30) . If the system 
(30) is examined in operator form, then the operator will act from L2 (where 
^ )?g - £ —dimensional Euclidean space, i,e., 


13 



Thus, we- may readily see that the operator which is conjugate to p"'*. 

will act from in and we -may show that the operator 


is the integral operator with the kernel 

i 




.6^ 

If we now take K=il as the operator K^gK s then equation (24) may he written 
in the form 


■ ^ *J r u \ <j / _ ' J' 

In its turn, equation (31) has the following analytical solution: 




• ■ i 




i \ 

where d = dg ) is found from the 


inversion of the matrix 


(31) 


(32) 


i. -tS) or 


(A+g) ^ (<?y j 


where 




a 


o:^ 


^ f 
/ 


i'5' 


Formula (32) has a very simp!^ _foim.^ ^ However, the analytical solution may 
also be obtained for the operator i ^ , We should also, note that there 

are methods making it possible to have a system of integral relationships and equa- 
tion (28). Thus, the integral relationships apparently are the finite form to 
which the systems of nonlinear parabolic equations are reduced. If we take a cer- 
tain more complex operator, as compared with that examined here, as the operator 
K = K, then in this case it is first advantageous tO' reduce the Fredholm operator 
to matrix form and then to use the algorithm (24) in the projection on the grid 


space. 



6. CONSERVATITO SINGLE STEP ALGORITHM Cll) - Cl2) , USING A KNOWLEDGE OP THE 
SYSTEM PREHISTORY. 


The iteration algorithm examined helow must be effective when is a 

bounded, poorly defined operator, When is a non-bounded operator in 

for example, the differentiation operator with respect to the evolutionary variable, 
the iteration process is directly used. We shall take its conservation to mean 
conservation within .an accuracy of the system approximation. However, with a 
knowledge of the system prehistory, which supplies certain a priori information, 
a single step -process is avoided. The quantities f and D are retained in the com- /18 
puter memory not only .in the n-1 layer, but also at certain previous layers up to 
n-S inclusively. It is natural .to assume that the older is the derivative of the 
function, the smoother it is. This is the basic a priori assumption here. Then, 

Instead of the process (13) - (14), we may regard the system' (11) - (12), where 
instead of equation (12) 



we shall solve the equation 


4 L(f 




The derivatives of all the functions used here may be calculated using the follow- 
ing formulas — •- . 

O'"'- 

■- - . ' 

where we set 0 = . It is apparent that the better known are the deriva- 

tives, the greater the amount of the. prehistory for solving the system, which is 
retained in the computer memory. We should note that if the dependence 
is strongly nonlinear, namely 

s' 0 

where , is a smooth function of t . » this casei)all the methods of di- 
rectly calculating for the differential operator "L give an increasing error, 

which disturbs the conservation. Let us examine the approach making it possible 
to retain a solution which is constant with respect to the error even in the case 
of strong nonlinearity. For this purpose, we turn to conjugate equations for which 



15 



the following is valid; 


Theorem 6 . The relative error of the diffusion coefficient is retained when 
changing to conjugate equations. This control of the basic system differs primari'- 

• A 

ly from the control of that examined above related to the addition of . The 
change to conjugate equations for the diffusion coefficient contains two basic ap- 
proaches, namely, the change to an implicit characteristic and then the change to 
the conjugate equation. If the operator , then 


or 



” "ft 

" “A+ ^ 
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Below, the indices a, ’3 ate omitted for purposes of simplicity. Then 

(33) 

Now let us multiply C33) by a certain function E Coagulator) 


(34) 

We require that the equation is satisfied. Then the equation (33) 

is equivalent to the following equation: 





(35) 


We should note that the diffusion coefficient D in (35) may be taken from the pre- 
ceding evolutionary layer, since the function E decreases and, consequently, its 
error does not increase. In turn, to calculate E we must use an implicit scheme. 
The accuracy of determining it will be , inhere . For this 

purpose, let us introduce the implicit characteristics (Fig, 1), 



Fig, 1 
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In conclusion, we should note that in connection with Fig. 1, the charac- /2Q 
teristics may he considered according to an implicit scheme, if their directions 
only depend on and simultaneously we use the step A^i. as the step At. 


7. ALGORITHMS FOR SOLVING THE SYSTEM OF 'EQUATIONS (7) - (8), USING THE 
APPROXIMATION OF HIGH ACCURACY,. 


As was shown above, the conservation requirement • Imposed on numerical al- 
gorithms for solving systems of nonlinear parabolic equations; is reduced to com- 
plex computational schemes with iteration control. However, in certain cases this 
requirement- may be reduced, if we solve the equation (11) very- precisely, for 
example, using the introduction of non-control difference grids or grids with 
"frequency increase," if this does not lead to a great increase in the amount of 
operations. This can be done if when solving equation (11) we use the method of 
variable directions. Briefly, we shall show how to retain the advantages of the 
method of variable directions for an increased frequency difference grid (Fig. 2) . 
These variable directions consist of reducing the numerical solution to satisfying 
the subsequent "trial run,” 


% 


«*■ 


Fig. 2 


For simplicity, we shall assume that has a diagonal form. Let us examine (11) 

■»/- . I'i; 

■Q-r 7 ^ ' 

Its difference analog in finite dimensional space, produced by the grid' in the 
domain Q, will be a linear system which within the framework of the variable direc- 
tion method may be written by the recursion relationship 






(ZVu,n)l“‘'= 


»wi 
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Generalizing the method of variable direction to a class of non-regular 
difference grids was jperformed in [33,- where it was, shoum how to separate the 
main operator into the operators ; and i for grids with 

partial frequency increase, in a two-dimensional domain. The' method introduced in 
[3] makes it possible to reduce- all the operator inversions -to satisfying the sub- 
sequent trial runs, with, a small increase in the number of equations. 


8. SCHEMATIC CLASSIFICATION OF NUMERICAL ALGORITHMS -FOR SOLVING SYSTEMS OF 
NONLINEAR PARABOLIC EQUATIONS, 

- No\ 7 let us briefly give certain practical results of' applying the approach . 
given above .to solving systems of equations with nonlinear diffusion. Since this 
problem is very cumbersome, each of the -algorithms proposed for the solution may 
not be used completely, but only partially. For, example, we shall' asstone that the 
dependence of the diffusion coefficients^ on' the solution is known. In 

this case, we should be able to expand >4^’ in Fourier series. However, another 
approach is advantageous .here, namely --i-in accordance with, the laws of conserva- 
..tion, -r_ irom the parabolic system C2).,. we--separate the retained quantities such- 
as energy, energy flux, etc. Then.*^^-' may be subjected to iteration, so that they 
satisfy the laws of conservation, .In conclusion, we would like to give- a certain 
schematic classification of. the numerical algorithms, which may be useful xirhen 
selecting a method for solving the systems (2) - C5) . 
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tig, 3, Practical diagram for selecting, the algorithm for 
the nxmerical solution of systems of nonlinear 
parabolic equations. 


In conclusion, we would like to express our appreciation to 
V. Gf Zolotukhin for assistance in the work and valuable comments, and we would 
also like to thank V. G, Popov for useful discussions, 
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